{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from popsycle import synthetic\n",
    "from popsycle import utils\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "In additional to the usual python/astropy things, you'll need to install Galaxia, PyPopStar, PopSyCLE.\n",
    "For PyPopStar, it goes faster once you've made isochrones.\n",
    "Same goes from PopSyCLE with the initial-final cluster mass ratio file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "# Galaxia\n",
    "\n",
    "This is where you generate the stellar model and set the survey direction and area.\n",
    "\n",
    "Galaxia requires a text file containing the parameters you want for the model,\n",
    "following the format given in the Galaxia example. That file can either be \n",
    "generated by hand, or with the `write_galaxia_params` function. \n",
    "\n",
    "To create the file `example_galaxia_params.txt`, run:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "synthetic.write_galaxia_params(output_root = 'example',\n",
    "                               longitude = 1.25,\n",
    "                               latitude = -2.65,\n",
    "                               area = 0.001)"
   ],
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "\n",
    "Then create the stellar model, \n",
    "by running `galaxia -r example_galaxia_params.txt` on the command line. \n",
    "PopSyCLE can also combine these two steps by executing:\n"
   ],
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "synthetic.run_galaxia(output_root = 'example',\n",
    "                      longitude = 1.25,\n",
    "                      latitude = -2.65,\n",
    "                      area = 0.001)"
   ],
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    " \n",
    "The results should look like this:\n"
   ],
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "\n",
    "CODEDATAPATH=/u/casey_lam/scratch/GalaxiaData/\n",
    "\n",
    "Reading Parameter file-             example_galaxia_params.txt\n",
    "___\n",
    "\n",
    "outputFile               example                 \n",
    "outputDir                ./                      \n",
    "photoSys                 UBV                     \n",
    "magcolorNames            V,B-V                   \n",
    "appMagLimits[0]          -1000                   \n",
    "appMagLimits[1]          1000                    \n",
    "absMagLimits[0]          -1000                   \n",
    "absMagLimits[1]          1000                    \n",
    "colorLimits[0]           -1000                   \n",
    "colorLimits[1]           1000                    \n",
    "geometryOption           1                       \n",
    "longitude                1.25                    \n",
    "latitude                 -2.65                   \n",
    "surveyArea               0.001                   \n",
    "fSample                  1                       \n",
    "popID                    -1                      \n",
    "warpFlareOn              1                       \n",
    "seed                     12                      \n",
    "r_max                    20                      \n",
    "starType                 0                       \n",
    "photoError               0                       \n",
    "___\n",
    "Reading tabulated values from file- /u/casey_lam/scratch/GalaxiaData/Model/vcirc.dat\n",
    "Using geometry:                     Patch at l ,b : (1.25 -2.65) d_theta=0.0178412\n",
    "Reading Isochrones from dir-        /u/casey_lam/scratch/GalaxiaData/Isochrones/padova/UBV\n",
    "Isochrone Grid Size:                (Age bins=182,Feh bins=34,Alpha bins=1)\n",
    "Time Isocrhone Reading              2.73553     \n",
    "Generating populations................\n",
    "___\n",
    "Thin Disk,ID=0:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_0_E1.ebf\n",
    "Time Tree generation/reading =      1.34589     \n",
    "Completed % <0..10..20..30..40..50..60..70..80..90..>\n",
    "Stars spawned =                     15          \n",
    "Time Spawning=                      0.519746    \n",
    "___\n",
    "Thin Disk,ID=1:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_1_E1.ebf\n",
    "Time Tree generation/reading =      0.894184    \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     785         \n",
    "Time Spawning=                      3.01123     \n",
    "___\n",
    "Thin Disk,ID=2:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_2_E1.ebf\n",
    "Time Tree generation/reading =      1.21969     \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     2304        \n",
    "Time Spawning=                      2.31239     \n",
    "___\n",
    "Thin Disk,ID=3:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_3_E1.ebf\n",
    "Time Tree generation/reading =      1.12221     \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     3299        \n",
    "Time Spawning=                      2.63636     \n",
    "___\n",
    "Thin Disk,ID=4:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_4_E1.ebf\n",
    "Time Tree generation/reading =      1.12869     \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     9040        \n",
    "Time Spawning=                      5.16435     \n",
    "___\n",
    "Thin Disk,ID=5:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_5_E1.ebf\n",
    "Time Tree generation/reading =      0.924357    \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     9983        \n",
    "Time Spawning=                      5.99943     \n",
    "___\n",
    "Thin Disk,ID=6:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_6_E1.ebf\n",
    "Time Tree generation/reading =      1.03473     \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     11927       \n",
    "Time Spawning=                      5.80737     \n",
    "___\n",
    "ThickDisk:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_7_E0.ebf\n",
    "Time Tree generation/reading =      1.00016     \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     6561        \n",
    "Time Spawning=                      1.83702     \n",
    "___\n",
    "Spheroid:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_8_E0.ebf\n",
    "Time Tree generation/reading =      0.320754    \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     589         \n",
    "Time Spawning=                      0.693446    \n",
    "___\n",
    "Bulge:\n",
    "Reading tree from file-             /u/casey_lam/scratch/GalaxiaData/BHTree-2.2/bhtree_with_wf/bhtree_9_E0.ebf\n",
    "Time Tree generation/reading =      1.0575      \n",
    "Completed % <0..9..19..29..39..49..59..69..79..89..99..>\n",
    "Stars spawned =                     620676      \n",
    "Time Spawning=                      28.1113     \n",
    "___\n",
    "Total stars written                 665179                  \n",
    "File written-                       .//example.ebf\n",
    "Calulating magnitudes................\n",
    "Reading Isochrones from dir-        /u/casey_lam/scratch/GalaxiaData/Isochrones/padova/UBV\n",
    "Isochrone Grid Size:                (Age bins=182,Feh bins=34,Alpha bins=1)\n",
    "Time Isocrhone Reading              2.67512     \n",
    "Calulating Extinction................\n",
    "Time for extinction calculation     0.876306    \n",
    "Total Time=                         73.89 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "# Population Synthesis\n",
    "\n",
    "The stellar model now needs compact objects, which we inject with PopSyCLE. \n",
    "This produces an .h5 file with both stars and compact objects, sorted by latitude and longitude.\n",
    "Here you can set the black hole and neutron star kick speeds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*********************** Starting popid 0\n",
      "Starting age bin  6.342803955078125\n",
      "Starting age bin  6.3\n",
      "Starting sub-bin  0\n",
      "Starting age bin  7.0\n",
      "Starting sub-bin  0\n",
      "Starting age bin  7.7\n",
      "Starting sub-bin  0\n",
      "Starting age bin  8.0\n",
      "Starting sub-bin  0\n",
      "*********************** Starting popid 1\n",
      "Starting age bin  8.176460266113281\n",
      "Starting sub-bin  0\n",
      "Starting age bin  8.384539353847504\n",
      "Starting sub-bin  0\n",
      "Starting age bin  8.592618441581726\n",
      "Starting sub-bin  0\n",
      "Found 2 stars out of mass range\n",
      "Starting age bin  8.800697529315947\n",
      "Starting sub-bin  0\n",
      "Found 18 stars out of mass range\n",
      "*********************** Starting popid 2\n",
      "Starting age bin  9.000540733337402\n",
      "Starting sub-bin  0\n",
      "Found 8 stars out of mass range\n",
      "Starting age bin  9.077965436220168\n",
      "Starting sub-bin  0\n",
      "Found 10 stars out of mass range\n",
      "Starting age bin  9.155390139102934\n",
      "Starting sub-bin  0\n",
      "Found 3 stars out of mass range\n",
      "Starting age bin  9.232814841985702\n",
      "Starting sub-bin  0\n",
      "Found 8 stars out of mass range\n",
      "*********************** Starting popid 3\n",
      "Starting age bin  9.301197052001953\n",
      "Starting sub-bin  0\n",
      "Found 10 stars out of mass range\n",
      "Starting age bin  9.347537861347199\n",
      "Starting sub-bin  0\n",
      "Found 13 stars out of mass range\n",
      "Starting age bin  9.393878670692445\n",
      "Starting sub-bin  0\n",
      "Found 21 stars out of mass range\n",
      "Starting age bin  9.440219480037689\n",
      "Starting sub-bin  0\n",
      "Found 15 stars out of mass range\n",
      "*********************** Starting popid 4\n",
      "Starting age bin  9.477176666259766\n",
      "Starting sub-bin  0\n",
      "Found 38 stars out of mass range\n",
      "Starting age bin  9.535049225091933\n",
      "Starting sub-bin  0\n",
      "Found 61 stars out of mass range\n",
      "Starting age bin  9.592921783924101\n",
      "Starting sub-bin  0\n",
      "Found 221 stars out of mass range\n",
      "Starting age bin  9.65079434275627\n",
      "Starting sub-bin  0\n",
      "Found 155 stars out of mass range\n",
      "*********************** Starting popid 5\n",
      "Starting age bin  9.698973655700684\n",
      "Starting sub-bin  0\n",
      "Found 99 stars out of mass range\n",
      "Starting age bin  9.73796232175827\n",
      "Starting sub-bin  0\n",
      "Found 125 stars out of mass range\n",
      "Starting age bin  9.776950987815855\n",
      "Starting sub-bin  0\n",
      "Found 125 stars out of mass range\n",
      "Starting age bin  9.815939653873443\n",
      "Starting sub-bin  0\n",
      "Found 38 stars out of mass range\n",
      "*********************** Starting popid 6\n",
      "Starting age bin  9.845118522644043\n",
      "Starting sub-bin  0\n",
      "Found 91 stars out of mass range\n",
      "Starting age bin  9.886338176012039\n",
      "Starting sub-bin  0\n",
      "Found 110 stars out of mass range\n",
      "Starting age bin  9.927557829380035\n",
      "Starting sub-bin  0\n",
      "Found 49 stars out of mass range\n",
      "Starting age bin  9.968777482748031\n",
      "Starting sub-bin  0\n",
      "Found 31 stars out of mass range\n",
      "*********************** Starting popid 7\n",
      "Starting age bin  9.94097840309143\n",
      "Starting sub-bin  0\n",
      "Found 158 stars out of mass range\n",
      "*********************** Starting popid 8\n",
      "Starting age bin  10.038501205444335\n",
      "Starting sub-bin  0\n",
      "Found 19 stars out of mass range\n",
      "*********************** Starting popid 9\n",
      "Starting age bin  9.9\n",
      "Starting sub-bin  0\n",
      "Found 17651 stars out of mass range\n",
      "Total run time is 52.840643 s\n",
      "******************** INFO **********************\n",
      "Total number of stars from Galaxia: 665179\n",
      "Total number of compact objects made: 43596\n",
      "Total number of things binned: 708775\n"
     ]
    }
   ],
   "source": [
    "synthetic.perform_pop_syn(ebf_file = 'example.ebf', \n",
    "                          output_root = 'example',\n",
    "                          iso_dir = '/u/casey/scratch/work/microlens/popsycle_test/isochrones/',\n",
    "                          bin_edges_number = None, \n",
    "                          BH_kick_speed_mean = 100, \n",
    "                          NS_kick_speed_mean = 350);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "# Finding microlensing events\n",
    "\n",
    "Here is when all microlensing events are found.\n",
    "You can set the survey duration, cadence, maximum impact parameter, blend radius.\n",
    "This also has the option of being parallelizable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on loop ll, bb =  0 0\n",
      "Total runtime: 33.141568 s\n"
     ]
    }
   ],
   "source": [
    "synthetic.calc_events(hdf5_file = 'example.h5', \n",
    "                      output_root2 = 'example', \n",
    "                      radius_cut = 2, \n",
    "                      obs_time = 1000, \n",
    "                      n_obs = 11, \n",
    "                      theta_frac = 2, \n",
    "                      blend_rad = 0.65, \n",
    "                      overwrite = False, \n",
    "                      n_proc = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "# Photometry\n",
    "\n",
    "This is the last thing, where you choose the photometric band for the observations and the reddening law.\n",
    "The final file of microlensing events is produced here as a .fits file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original candidate events:  117\n",
      "Candidate events in survey window:  99\n",
      "Total runtime: 0.580987 s\n"
     ]
    }
   ],
   "source": [
    "synthetic.refine_events(input_root = 'example', \n",
    "                        filter_name = 'I',\n",
    "                        photometric_system = 'ubv',\n",
    "                        red_law = 'Damineli16', \n",
    "                        overwrite = False, \n",
    "                        output_file = 'default')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "# Things you can do with PopSyCLE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "First, read in the table, and print out a list of all the possible keys:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['zams_mass_L', 'rem_id_L', 'mass_L', 'px_L', 'py_L', 'pz_L', 'vx_L', 'vy_L', 'vz_L', 'rad_L', 'glat_L', 'glon_L', 'vr_L', 'mu_b_L', 'mu_lcosb_L', 'age_L', 'popid_L', 'ubv_k_L', 'ubv_i_L', 'exbv_L', 'obj_id_L', 'ubv_j_L', 'ubv_u_L', 'ubv_r_L', 'ubv_b_L', 'ubv_h_L', 'ubv_v_L', 'zams_mass_S', 'rem_id_S', 'mass_S', 'px_S', 'py_S', 'pz_S', 'vx_S', 'vy_S', 'vz_S', 'rad_S', 'glat_S', 'glon_S', 'vr_S', 'mu_b_S', 'mu_lcosb_S', 'age_S', 'popid_S', 'ubv_k_S', 'ubv_i_S', 'exbv_S', 'obj_id_S', 'ubv_j_S', 'ubv_u_S', 'ubv_r_S', 'ubv_b_S', 'ubv_h_S', 'ubv_v_S', 'theta_E', 'u0', 'mu_rel', 't0', 't_E', 'ubv_i_app_S', 'ubv_i_app_L', 'cent_glon_i_N', 'cent_glat_i_N', 'ubv_i_app_N', 'ubv_i_app_LSN', 'delta_m_i', 'f_blend_i', 'pi_rel', 'pi_E']\n"
     ]
    }
   ],
   "source": [
    "t = Table.read('example_refined_events_i_Damineli16.fits')\n",
    "print(t.colnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "## Figure out how many events due to BH, NS, WD, star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of BHs: 2\n",
      "Number of NSs: 1\n",
      "Number of WDs: 15\n",
      "Number of stars: 81\n"
     ]
    }
   ],
   "source": [
    "bh_idx = np.where(t['rem_id_L'] == 103)[0]\n",
    "ns_idx = np.where(t['rem_id_L'] == 102)[0]\n",
    "wd_idx = np.where(t['rem_id_L'] == 101)[0]\n",
    "st_idx = np.where(t['rem_id_L'] == 0)[0]\n",
    "print('Number of BHs: ' + str(len(bh_idx)))\n",
    "print('Number of NSs: ' + str(len(ns_idx)))\n",
    "print('Number of WDs: ' + str(len(wd_idx)))\n",
    "print('Number of stars: ' + str(len(st_idx)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "## Histogram of tE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'$t_E$ (days)')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAE0CAYAAAAhaTThAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEvxJREFUeJzt3X2MZXV9x/H3p1JFaUPBRY0P64gjID7UyjSlmLa0ii3a9SESGluJ2MpijBprql0tTZrYB5oYbGpFpaZFrY2KWHWLiVWgC6JI2GrV9QEU1yekLoJQRCzKt3+cM+043LvM7JyZe+/83q/k5s79nXN+53v3t5nPnOdUFZKkdv3UpAuQJE2WQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktS4gyZdwFpt2bKl5ubmJl2GJE3c7t27b6yqI1a73MwHwdzcHFdfffWky5CkiUvytQNZzl1DktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMbN3AVlSbYB2xY/z8/PT7AaqTO346LB+9x79tMH71MaZeaCoKp2AjsXPy8sLJwxwXIkaea5a0iSGmcQSFLjDAJJapxBIEmNMwgkqXEzd9aQtFbrcaqnNMvcIpCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXFeRyBNKW9trY3iFoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4yYSBElek+RLSe5K8qxJ1CBJ6kxqi+Bi4GnAZRNavySpd9AkVlpVnwRIMonVS5KW8BiBJDVuIlsE0krN7bho0iVsKuvx77n37KcP3qc2llsEktQ4g0CSGjep00fPSvJN4JeBtyb5ZpIHTaIWSWrdioIgyUOTvCHJJ5LcnqSSzI2Z92FJ3pvkliS3Jnlfkq1L56mqP6+qh1bVfapqS//zDWv/OpKk1VrpFsE8cCpwM3D5uJmS3A+4BDgGeD5wGvAo4NIkh6ytVEnSeljpWUOXVdUDAZK8EHjqmPnOAI4Ejq6qL/fzfwa4FjgTOGe1BSbZDmwfN33r1q3jJkmSVmBFQVBVd62wv2cAVy6GQL/sV5NcATyTAwiCqjoPOG/c9IWFhVptn5Kk/zf0weLHAJ8b0b4HOHbgdUmSBjB0EBxOdxxhuZuAwwZelyRpAOtx+uioXTXeVEiSptTQQXAz3VbBcocxektBkjRhQwfBHrrjBMsdC3x+4HVJkgYw9E3nPgi8LsmRVXUdQH/h2ZOAHUOsIMk2YNvi5/n5+SG6laRmrTgIkpzS/3hc/35ykn3Avqra1bf9PfAS4ANJzqI7XvBa4BvAW4YouKp2AjsXPy8sLJwxRL+S1KrVbBFcsOzzuf37LuBEgKr6fpLfAF4PvIPuIPHFwMur6ra1lSpJWg8rDoKqWtGZP1X1deA5B1yRJGlDeRtqSWqcQSBJjTMIJKlxBoEkNW7mHl7vdQSSNKyZCwKvI5CkYblrSJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxs3cdQReUDaMuR0XTboESVNi5oLAC8okaVjuGpKkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMbN3JXF3mJCkoY1c0HgLSYkaVjuGpKkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUuJm7xYT3GpKmy9yOiwbvc+/ZTx+8T403c0HgvYYkaVjuGpKkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWrczN191NtQS9KwZi4IvA21JA3LXUOS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkho3c4+q9JnFkjSsmQsCn1ksScNy15AkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNe6gSRewWkm2AdsWP8/Pz0+wGkktm9tx0eB97j376YP3eU9mLgiqaiewc/HzwsLCGRMsR5JmnruGJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxk0kCJI8MsnHklyT5FNJFiZRhyRpclsEbwbOr6qjgFcB70ySCdUiSU3b8CBIcgRwPPA2gKr6SD/puI2uRZI0mS2CrcD1VXXnkrav9e2SpA02LQeL3S0kSRMyiSD4OvDgJD+9pO3hfbskaYNteBBU1T7gKuB0gCQn0W0R7N7oWiRJcNCE1vsi4G1JXgncDvxeVdWEapGkpq1oiyDJQ5O8IcknktyepJLMjZn3YUnem+SWJLcmeV+SnzgQXFXXVtUJVXVUVT2hqq5a+1eRJB2Ile4amgdOBW4GLh83U5L7AZcAxwDPB04DHgVcmuSQtZUqSVoPK901dFlVPRAgyQuBp46Z7wzgSODoqvpyP/9ngGuBM4FzVltgku3A9nHTt271rFNps5nbcdGkS2jKioKgqu5aYX/PAK5cDIF+2a8muQJ4JgcQBFV1HnDeuOkLCwseW5CkNRj6rKHHAJ8b0b4HOHbgdUmSBjB0EBxOdxxhuZuAwwZelyRpAOtxHcGoXTVeOSxJU2roILiZbqtgucMYvaUgSZqwoYNgD91xguWOBT4/8LokSQMY+sriDwKvS3JkVV0H0F949iRgxxArSLIN2Lb4eX5+fohuJalZKw6CJKf0Py4+N+DkJPuAfVW1q2/7e+AlwAeSnEV3vOC1wDeAtwxRcFXtBHYufl5YWDhjiH4lqVWr2SK4YNnnc/v3XcCJAFX1/SS/AbweeAfdQeKLgZdX1W1rK1WStB4y6/d667dKvtZ/PBS4ZcRso9q3ADeuY2mrNa72SfS52uVWMv89zbO/6Y7rMP05rsOZ1nF9eFUdseoeqmrTvIDzVtoOXD3peldS+yT6XO1yK5n/nubZ33THdZj+HNfpGIdJj+uo17Q8oUySNCEGgSQ1brMFwc5Vtk+T9ajxQPtc7XIrmf+e5tnfdMd1mP4c1+HM8rjezcwfLD5QSa6uqoVJ16FhOa6bk+O6vjbbFoEkaZUMAklqXMtBMPZhN5ppjuvm5Liuo2aPEUiSOi1vEUiSMAgkqXkGwRhJHpnkY0muSfKpJJ66tgkkeU2SLyW5K8mzJl2P1i7JwUnen+QLST6d5MNJjpx0XbPEIBjvzcD5VXUU8CrgnUl85Obsuxh4GnDZpAvRoN5UVY+uqifQXWD11kkXNEsMghGSHAEcD7wNoKo+0k86buxCmglV9cmq+sqk69BwquqOqvrwkqYrAbcIVsEgGG0rcH1V3bmk7Wt9u6Tp9lLgA5MuYpYM/ajKzczdQtKUS/Jq4CjgyZOuZZa4RTDa14EHJ/npJW0P79slTaEkfwQ8Bzi5qm6fdD2zxCAYoar2AVcBpwMkOYlui2D3BMuSNEaSVwDPBU6qqu9Nup5Z45XFYyR5FN3B4i3A7cD2qrpqslVprZKcBbwIOAL4b+AOYKGqbphoYTpgSR4KfAO4jm5MAX7k3UpXbtMEQf+f4Y+BBeDngfsCj6iqvSPmfRjwemDxL/2PAi+vKnf9TBnHdXNyXKfLZto1NA+cCtwMXD5upiT3Ay4BjgGeD5wGPAq4NMkhG1CnVsdx3Zwc1ymymc4auqyqHgiQ5IXAU8fMdwbdOcZHV9WX+/k/A1wLnAmcswG1auUc183JcZ0im2aLoKruWuGszwCuXPxP1S/7VeAK4JnrUZsOnOO6OTmu02XTBMEqPAb43Ij2PcCxG1yLhuO4bk6O6wZoMQgOp9svudxNwGEbXIuG47huTo7rBmgxCABGnSrllcOzz3HdnBzXddZiENxM91fGcocx+i8PzQbHdXNyXDdAi0Gwh26/43LHAp/f4Fo0HMd1c3JcN0CLQfBB4PilD65IMgc8qZ+m2eS4bk6O6wbYNFcWAyQ5pf/xyXS3EXgxsA/YV1W7+nkOAf4T+AFwFt3+x9cCPws8vqpu2+i6tX+O6+bkuE6PzRYE477Mrqo6ccl8W/nJS9Yvprtkfe9616jVc1w3J8d1emyqIJAkrV6LxwgkSUsYBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQQSkOTtSb6zmgeiJ9mbZO86lnXAkhyXpJL8waRr0fQzCDQTkryi/8X2u+vQ9wLwPODsqvr+0P1PQlXtBt4P/HmSn5l0PZpuBoFmxRP7993r0PdfArcCb1qHvifpr4AHAS+bdCGabgaBZsVxwG3AtUN2muQo4CnAe6rqB0P2PWlVdRXwReDMJPeadD2aXgaBplqSv+5vV3wM8DPAj/tdRJXkeQOs4vfpbm387jHrT5KXJNmT5I4k30ryd0kO3U/Npye5MMl1SX6Q5NYkVyyvN8kx/fe4ZD99fTbJnUketKTtGUkuTvLtJD9Mcn2SXUlePKKLdwFb6cJOGumgSRcg3YPdwNuA5wMfBz6yZNquAfp/CvBj4Mox0/+GbtfKt4HzgDuBZwK/BNwb+J8Ry7yJ7jGKl/XL3R94GvCOJEdX1Z8CVNUXk1wK/HqSo6rqmqWdJDkBeCxwYVXd0LdtB94C3ADsBG4EHgA8HngBcO6yWq7o308CPnxP/xhqVFX58jXVL2A73ZOptg/c7yHAj4DPjpl+Qr/eLwOHL2k/GPhEP23viOUeOaLt3nQPVLkTeMiS9lP6fl43Ypnz+2knLWnbDfwQeMCI+beMaDu07+OqSY+jr+l9uWtIs2DxQPF/DNzvQ4B70f3VPsoL+ve/qKqbFhur6g7g1eM6raqvjGj7H+CNdFvhT14y6f3A9cDpSe6z2Jjk54BTga8AH13W3Y/oAmX5Om4c0XYLcAfd7iFpJINAs+CJdL/4Prt8QpJLlhwzWPq6aAX93r9/v3k/64XRu6Aup/uFfDdJtiZ5Y5IvJrl9sSbgwn6WhyzOW1U/At7a1/KcJd2cBtwXOK+qlj5G8J3A/YA9SV6f5FlJjhj/FQG4CdhyD/OoYR4j0FRLchDwOODzVfXDEbP8AvAnwD8sa799Bd0vniV08JjpiweE/2v5hKr6cZLvjqj3SOAq4DC6sPg34Ba64xBzdMc67rNssfOA1wBnAv/ct22nO/7wj8vWe06SG+ke9P4y4OVAJdkFvLKqrh7xPe675LtKd2MQaNodS/eL+m67hZI8Evg5uoed33AAfX+nf7//mOm39O8PBK5btu579ct9a9kyr+jbX1BV5y9b5rl0QfATqupbSXYCz07yaLoQeSzw7qraN2L+twNv73cfnQA8m+7spw8neXRVLX4vkvwU3b/RV8d8R8ldQ5p6T+jfPzVi2nF0f2mPmrYS3wb2AUePmb4YPr82YtqvMPoPqfn+/cIR00b1s2jxbJ/t/Qu6s4PGqqrvVdWHquoMugPLh/d1LXU03emxn95fX2qbQaBpt/jX+q0jph1Hd7D3O0luW/K6YCUd9/veLwO2JJkfMcv5/fufJDl8sTHJwXRX7Y6yt38/cWljkt8EXrifci4GrqHbYjgVuKaqLl0+U5Lf6neXLfeA/n35LrHj+/e79SUtcteQpt3iLSX+Isljge8De6rqAroguBDYsWyZW1i5C+kO0v4m3Wmi/6eqrkjyBuClwOeSvJf/v47gZkafbXQu3dlGFyS5kG7X0WOB3wLeA/zOqCKqqpK8GTinbxq3NfAu4I4kH6MLndBtBfwi3b/V8jOMnkq31fSBMf1JXkfga/pfwEuAL9GdBll0p3MCfBf4wzX2fW+6i7M+OWZ6+vV/ge78/evpTgM9lO4X8d4Ry5wAXEIXFv8NfAx4Ft1WQgF/NmZdh9H90r4DuP+YeV4E/AvdMYvb6c4I+hTwKuBnl817KN1B4vdPegx9TfcrVUvPTJNmQ5JH0P0y/PWq+vc19vVquhvPPbGqDvR4w5olOZFuF84/VdVpA/T3UuBvgV+tqsvX2p82L4NAMynJKcAFwKOB7y2bfGN15+evtK+D6bY4PlNV24arcnWSfAg4GTi+qj65xr7uS3cx2ser6pQh6tPm5TECzarj+vcvLGsvurNnlofDWFV1R5LT6O75c0ht4DMJkjwO+G2673My8K9rDYHeHN31CecP0Jc2ObcIpAlKcjrdRWO30t0U7sU14lYR0noyCCSpcV5HIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWrc/wLP9CjVo4RhAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(t['t_E'], bins = np.logspace(0, 2.5, 15))\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('$t_E$ (days)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "## piE vs tE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'$\\\\pi_E$')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAE0CAYAAACB9lwPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X30HFWd5/HPlxAg4G4CgnoMZgJGw5NC4LdnWLMPGhVGmWBGGOa4MxxklwSPq3tmmcmZoOwOjDhmDopzxh1YGXfX1XEdhLjB3+I52ZVgwIzgJoYHoyDKoz90BUnCQBLMw3f/qOqkf52q7uruqrr18H6d0+eX7q6uut1dqW/fe7/3XnN3AQAQ0mGhCwAAAMEIABAcwQgAEBzBCAAQHMEIABAcwQgAEBzBCAAQHMEIABAcwQgAEBzBCAAQHMEIABAcwQgAEBzBCAAQHMEIABDc4aELUHfHH3+8z58/P3QxACC4zZs3P+/uJ4zyWoLRmObPn69NmzaFLgYABGdmT436WprpAADBEYwAAMERjAAAwRGMAADBEYwAAMERjAAAwRGMAADBMc4IqKi1W6Z0w7pH9ez2XXr9nFlaef5CLVs0N3SxgEIQjIZkZkslLe3cX7BgQcDSoKnWbpnS1V9/WLv27JMkTW3fpau//rAkEZDQSDTTDcndJ919Rec2e/bs0EVCA92w7tEDgahj1559umHdo4FKBBSLYARU0LPbdw31OFB3NNPVGH0KzfX6ObM0lRB4Xj9nVoDSAMWjZlRTnT6Fqe275DrYp7B2y1TooiEHK89fqFkzZ0x7bNbMGVp5/sJAJQKKRc2opvr1KVA7qr/Od0jNF1nVvaWEYFRT9Ck037JFc2t1MUE4Tci+pJmuptL6DuhTANqnCdmXBKOaok8BQEcTWkpopqsp+hRGM0y7et3b4NEeTci+JBjVGH0KwxmmXb0JbfBoj5XnL5x2vkr1aymhmQ6tMUy7ehPa4NEeyxbN1afe/xbNnTNLJmnunFn61PvfUqsfTtSM0BrDtKs3oQ0e7VL3lhJqRmiNYTIQyVYEykUwqpG1W6a0ePV6nbTqTi1evb7Rsy0U8V6HyUAkWxEoF810NdGmDvWi3uswGYhkKwLlMncPXYZaSVjPaPljjz1W+HEXr16fmLo5d84sbVy1pNBjl53iHPK9Vg3p5agTM9vs7hOjvJaa0ZDcfVLSZOf+xMTE8jKOG6pDPUSNjOSBSJtqwwB9RjURqkM9RIozyQMR0svRJgSjmgjVoR6ilkLyQIQaItqEYFQToQa1hailNGEAXx6oIaJNSGAY08TEhG/atCl0MQYatSO8t99CimopbQwOZUv67CVpzqyZuvbC0/n8UTkkMKCvcTrCSXEOp/MZXze5Vdt27jnw+PZde0hkQONQMxpTHWpGpErXWxO/P1LWm4maEfqiI7zemvb9kbKOJASjFmjCWidtluf3l0eNZNx99EtZJxi1F8GoBZqw1kmb5fX9jVMj6QSgqe27ZJI6jfuj1GqaVtNDPkjtboEsqdJtmoS1bvJKdR91EG0niHVqZ729zMMOxCVlHUmoGbVEv7VOaMOvvjzWqhm1RpIUxIbdRzdq6khCzQhMO9MSo9ZIsgSaYWo1DGpGEmpGoA2/JUatkaQlUAyzj151X5UU+aNmhMq14dN/VYxRayRJcwVa/JdaDfJCzQiVasNvav9VVQZ5jlIjYRYOlIFghEpdbJo4BqUJAZZmNRSNYARJ1bnYNLH/qokBFsgbwQiV0sTZIsYJsFVp3gOKRgIDKqWJC+uNmiDSPdjUdbB5j4QONBE1IwTT71f/oNpAv9dWrTYxaoIIzXtoE4IRghjUqd/vYtvvtZIqlywwaoJIE/vPgDSsZzQkM1sqaWnn/oIFC5Y/9thjAUtUT8Ou0dNd2znMTPsSztu5cbNX0n7nzJqpY448PFMwKKtmNeg4TVzHCM3GekYlcvdJSZOd+xMTE8sDFieTKjRb9ZYhbUR/0q/+3ppQUiBKe23H9l17tH1XtFpqv9pSWWnYWY5TpfFfdVSF8x7ZEYwablCTVlk1gN4ydC9D0C2pUz/LRJ3dr+03dU1HWt9LWf00WY5T9vivJl28mzC2q20IRg2XdtG79htb9cre/aX8Z00qg0uHBKSZM0wvv7JXJ626c9rFMEsfSXeNobc2kSZpv2X102Q9Tlnjv6rwoyVPJH/UD8Go4dIuep0mq25F/WdNK4Mr6v94dvsuzTl6pl7avTexKS2tWW+Gmfa7J14guy+eO3+9V9t2Hvp+k2phacc6zOyQIDmOqo2nqsKPljyR/FE/jDNquGEvbkX8Z00rQ6cj/onVF+joIw7Xnv3TG+46wTFt7NFnLjlTT6y+QBtXLZl2YVy2aO6B/W5ctUR/uvT0zGOXko4lRf1UeY71qdp4qn4/Wuq4vEjVJv/FYASjhku76B179MzE7Yv4z5rlwtvvl+y4698M8/rebWeYHbJNHhfjqq3pU4UfLXmqWrDHYDTTNVxaJ7h0aN9KUf9Zs3TED2q2GqfvZNiO+e5jnbTqzsRt8rgYV2U+QCk9c++omYdlbuKskipN/otsCEYt0O+iV9Z/1kEX3pXnL9TK2x6c1lQ38zAbOziOm1VVtb6dolThR0veqhTsMRjBqMUq95+1t0Xs0BayoY2bVdWmsT5V+NGC9iIYoRJuWPeo9uybnsCwZ5+Pnd03blYVzT0V/NGCRiIYoRKKSsXNo5mtbRfjJg1+RX0QjBBE7wVvztEzC+koH6eZrY0X5ZAzF7Tx88ZBBCOMZJwLR9IFb+ZhppkzbFpTXVrQGObYozaztXU6mVAzF7T188ZBBCMMbZQLx6BZt/fs90wza49y7FGa2Yq+KFe1FhBq5gKm7wHBCEPLeuHoXHB7J0ZNm3V7x649euBPz8vl2OMq8qJc5VpAvz62IgMo0/eAGRgwtCwXju4ls6XkGbp7ZekfKuuiVeR0Mv0CamhpMxe845QTCl0Cnel7QDDC0LJcOLIu+9CRNamgrItWkdPJVLkWkDZN0d2PPFdoAGX6HtBMh6FlyVDLcmHtN+v2OMfOQ5Hji6o+q0NSH9u/v/WBxG3zCqCM5wLBCEMbZ665jlkzZ4w0MWiZF63ei/LaLVNavHr92Met46wOZQTQto3nwnTmKZ3JyGZiYsI3bdoUuhiV09tJLx1cTG9uDX/1Jr2fUQNqZ391qgXk/f7RTGa22d0nRnktNSMUomnNLnln8dWtFtC07xPVQ81oTKPUjOr2q7itur+ntP8lJumJ1ReUWSygsqgZ1UiVx5jgoKRmqSRFJx004YdLE94Dikdqd8mqPMYEB2VJTS866aB7rFYRY3vK0IT3gHIQjEpW5TEmOKjf91HWMuFN+OHShPeActBMV7KqjzFpqmGbitK+p87YqDI04YdL1vdAUx6oGZWMkeblG6WpKOl7kqJ59cpqbmrCFDlZ3gNNeZAIRqVLm26FX4HFGaWpqPd7mmGHroFedHNTE364ZHkPNOVBanEznZl9TNJlkt4k6f3uvrasY9dtjEndjdrc1f09nbTqzpH2MY4mjO3J8h6a0ByJ8bU2GEm6S9Ktkv5L6IKgWHn004Xq62vCD5dB74F+VEgtbqZz9/vd/aehy4Hi5dHcFbLJrDMn3kmr7tTi1esb15fShOZIjK/NNaPaIvNoOHk0dy1bNFebnnpBX73/Ge1z1wwzXXRO8bWWa9Y+rK/c9/SBGSCaOEi6Cc2RGF/jgpGZfV/SvJSnF7n7M2WWJ2/M4DCacZu71m6Z0prNUwdWqd3nrjWbpzTxG8cV9rmv3TI1LRB1NHE57iY0R2I8jQtG7n526DIUqaxlt5tu2NpliM/9hnWPps6JR+c+mqZxwajpyDwaXm/geccpJ2jN5qmhapchPvd++6ZzH03T2gQGM7vGzH4m6Z9K+oKZ/czMXhe6XIM0YSBkmZIGVH7lvqeHHtcS4nNP27dJdO6jcYIHIzM70cw+Z2bfNbOdZuZmNj9l2zeY2e1mtsPMXjSzr5tZWv9QX+5+vbuf6O5Huvvx8b9/Mc57KQOZR8NJal4bpekrr899mMy4pGOapN8/d16lm2Sbnv3XVKG/tyo00y2QdImkzZLulXRe0kZmdrSk9ZJeUTRY1SVdL+luM3uru79cTnHD6lyErpvcqm0790iSjjw8+G+KyhqmGa1fLSePjK9hk0/KzDLLK0OTBJt6qsL3VoVgdI+7v1aSzOwKpQQjScslnSxpobv/JN7+IUmPSbpS0o1FFM7MVkhakfb8vHkjVczGtnvP/gP/3r5rD//hU6QNqOwsgd6RpZYzbsbXKEkQg46ZRxDJ80JEgk09VeF7C/6T2t33D95KknShpPs6gSh+7ROSNkp6XxFli49xi7tPpN1OOOGEog6dirm8sktrXvv9c+eVPj9g3kkQgyYYzdrskuf5RIJNPVXhe6tCzSir0yXdkfD4Vkm/W3JZgqrCiVMXVRpQmfe0N4OCSNbaTp7nE1P71FMVvrfgNaMhHCdpW8LjL0g6tuSyBEVG3XCWLZqrjauW6InVF2jjqiXBmovyTj7pF0SGqe3keT6RYFNPVfje6hSMpOREqEPn9m+4Kpw4TVFEBlHaPvNePqRfEBmmtpPn+cQSKfVUhe+tTs102xTVjnodq+QaU2NVqempzorIIFq7ZUorb39Qe/b5gX2uvP3BA/vMc9qblecvnFZ+6WAQuWHdo5mbXfI+n5jap55Cf291CkZbFfUb9TpN0g/LKoSZLZW0tHN/wYIFZR16mtAnThMUkUF03eTWA4GoY88+13WTW4dOAx8UHAYFkbRAlYTzCaHVKRh9Q9Knzexkd39ckuLBsYslrSqrEO4+KWmyc39iYmJ5WcdGvopIBOmM/cr6eJJhamxpQYTaM+qmEsHIzC6O/3lO/Pc9ZvacpOfcfUP82N9I+oikO8zsGkX9R5+Q9Iykz5dZXjRD6AyitNpPXjU2ajuok6okMNwW3z4U378pvn9dZ4N4hoUlkn4s6cuSviLpCUlL3P2lUkuLRigiEWTOrJmJj5tpWnJEvzFCpO6jjSpRM3L3TBlx7v60pIsKLk6jsTDfQUU0ZV174elaeduD2rN/er+Ru6Y1tfWr/YSusQ3COYQiVCIYoRxVmH+qavJuyurs64++9uCBhfg6upva+tV+Pvt7Zw2VfFAmziEUpSrNdCgB0wjlL2lM0bJFc7Xfk+cGn9q+S4tXr0+dOfz1c2ZVYsxHGs4hFIWa0ZCqkto9Cvoi8tWvltBvgtakx6XptZ+qJh9wDqEoQ9eMzGzSzF5VRGHqwN0n3X1F5zZ79uzQRcqMaYTy1a+WkLYWUVqNqEq1n344h1CUUZrp3ivp6M4dM7vVzF7ddf8wM/vHeRQO+WrjNEJFLhjWr5aQ1NSWFohMCjpn3jDaeA6hHKM00/Vmvr1X0mxJv4rvnyBpasR9o0BtGwhZdGf7oKy33qa2xavXVzpLLouQ5xBZfM1WVMAgMaKiqtoXUYSiFwzrNzdcHttXVYhziCy+5isqGKW1SAClydLZPs6v7WFrCW2rmeapCiuRolijBqPLzWyDpAfi+wQfVM6gZrQ8fm0PW0vIsj3NUYcii6/5RglG35b0J5I+JWlvvI+/MLONkr4v6Ze5la6C6pza3TaDmsWK/rWdFFQ6x00LNDRHJav6rBQYn3nK4LyBLzQ7WdHEpp3bIh1cb8jdfUbaa5tkYmLCN23aFLoYrZO19tBvu5NW3Zm6WuMTqy8Yu3y9gXDmYSaZpi0xMWvmjGkp3WlJDnPnzNLGVUvGKlOdJX2evZ8dwjOzze4+McprB9aMzOyt7v5Q7+PxMg6PK5rQtLPtfEkTks4epTBAFnkssSAV+2s7qdbVO1+ddGhNjOaoZPS3NV+WZrotZvYf3f2TgzZ09yclPSnp9jHLBaTKq3mtyOy2YYJH97Y0R6VrUyZoG2VJwTZJ55nZ5Wa2yMwOmSPfzOaa2VX5Fw84VF61hyLngBsmeHRvm2VQaZEDeYFQsiYw/LP4Jkl7zewRRZl0D0h6WNIZkv5M0o25lxC1UGYGWJ61h6J+bSfVutL6jLoDzbJFc7XpqRf01fuf0T53zTDTReccLCMJDmiqrINTb1W0tPfXFDXDnSHpUkmfkbQu/ruzgPKhBvotFFeEOkxJk1TruuF3z9QNF5/Ztya2dsuU1myeOrD8xD53rdk8deCzZNZsNFXWmtEj7n5D546ZHaMoe+5sRYHpOEWBCi1U9oDEunRmp9W6+pVz0GdJggOaaqRBr/ES4N+Jb2i5EBfIpnZmD/osSXBAUzGZ6ZAY9HooLpD5GfRZ1ml+O2aSqJfQ31fWPqOLzezPzOz9ZvbGQktUcXVez6godejDqYtBn2WVV4HtVnY/IsZThe8ra83ojPjmkmRmL0l6SNKDijLqHpT0sLvvLqKQqLa69OHUQVo2nRTNzhDy8x3mlzMTm9ZLFb6vrMHoJknfknRW121xfOvkqe6TdETeBUQ9NLUPp2xJ2XS3/t9ndOv3njkwg0OIdO5hU8pJtKiXKnxfWZrpflvSV919rbtf6+7L3H2+ogy6JZL+SNLfStpaXDGBdkicRmifHzKVUNnp3MOmlLM8eb1U4fsaGIzc/ZvuvjHh8e3u/m13/0t3v8zdFxVTRKBaQixlPu624xr2lzP9iPVShe+LbDpgCKGWMk/btizDZkzSj1gvVfi+Rl5CAhGWkGiXopd4SFx6YoZJPn3W77KXT2AJB2RR6BISQNv0yxoruqM37Rdq0mNlBoEq/HJGsxGMgC6DmuHKGOA7yjRCZSBjEkXKOugVMTNbama3dG47duwIXSTkaFDWWBU6eoEmomY0JHeflDTZuT8xMbE8YHGQs0HNcDRXAcUgGAFdsjTD0VwF5I9mOqBL3ZvhWAUWdUXNCOhS52Y4VoFFnRGMgB51bYarwmSXwKhopgMaogqTXQKjomZUQaEXuUI9scgh6oyaUcVUYZEr1FPdky/QbgSjihl2qn6goy6rwAJJaKarGNr9MY66Jl8A1IwqpgqLXAFA2QhGFUO7P4A2opluSGa2VNLSzv0FCxbkuv86D7oEgFGxuN6YWFwPACLjLK5HMx0AIDiCEQAgOIIRACA4EhgA1BLTZjULwQhA7bBcRvPQTAegdpg2q3kIRgBqh2mzmodgBKB2mDareQhGAGqHabOahwQGALXDtFnNQzACUEssl9EsNNMBAIIjGAEAgiMYAQCCo89oSEWvZwQAbUQwGpK7T0qa7NyfmJhYHrA4ANAINNMBAIKjZgQABWJ28WwIRgBQEGYXz45mOgAoCLOLZ0cwAoCCMLt4dgQjACgIs4tnRzACgIIwu3h2JDDUBBk5QP0wu3h2BKMaICMHqC9mF8+GZroaICMHQNMRjGqAjBwATUcwqgEycgA0HcGoBsjIAdB0JDDUABk5AJqOYFQTZOQAaDKa6QAAwRGMAADBEYwAAMERjAAAwRGMAADBEYwAAMGR2j0kM1sqaWnn/oIFCwKWBgCagWA0JHeflDTZuT8xMbE8YHEAoBFopgMABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAERzACAATXymBkZkeZ2Voz+5GZPWBm68zs5NDlAoC2amUwit3s7qe6+1mSJiV9IXSBAKCtWhmM3H23u6/reug+SdSMACCQVgajBB+VdEfoQgBAWx0eugB5M7PvS5qX8vQid3+mZ/urJb1Z0juLLhsAIFnjgpG7n511WzP7Y0kXSXqXu+8srlQAgH4aF4yyMrOrJH1AUSDaHro8ANBmrQxGZnaipM9IelzS3WYmSXvdfSJowQCgpYInMJjZiWb2OTP7rpntNDM3s/kp277BzG43sx1m9qKZfd3M0vqHUrn7z9zd3P2N7n5WfCMQAUAgwYORpAWSLpG0TdK9aRuZ2dGS1ks6RdJlki6V9CZFNZtjSignAKAgVWimu8fdXytJZnaFpPNStluuaCzQQnf/Sbz9Q5Iek3SlpBuLKJyZrZC0Iu35efOGrpgBAHoED0buvj/jphdKuq8TiOLXPmFmGyW9TwUFI3e/RdItac9PTEx4EccFgDapQjNdVqdL+kHC41slnVZyWQAAOapTMDpOUb9SrxckHVtyWQAAOapTMJKkpCYxK70UAIBcBe8zGsI2RbWjXscqucZUis2bNz9vZk+NsYvZknaU/PphXpN120HbDXr+eEnPZyxTVY37XVbhmJyPEc7H0bxp5Fe6e2Vukq5QVPuZn/DceknfSXj825I2lFjGpYoSGjq3pWPu75ayXz/Ma7JuO2i7DM9vKus7LPDcGOu7rMIxOR8PPM/5WPIx61Qz+oakT5vZye7+uCTFg2MXS1pVViHcfVLR+kcAgJxUIhiZ2cXxP8+J/77HzJ6T9Jy7b4gf+xtJH5F0h5ldo6gG9QlJz0j6fJnlBQDkqxLBSNJtPfdviv9ukPR2SXL3l81siaTPSvqyosSFuyT9obu/VFI5izBuLWuU1w/zmqzbDtquDbXJEO8x72NyPjZHrc5Hi9v5gKDMbJMzPyAqgvOxfHVL7QYANBDBCAAQHMEIVZE6/x8QAOdjyegzAgAER80IABAcwQgAEBzBCJVnZh8zs0fNbL+ZLQtdHrSXmR1lZmvN7Edm9oCZrTOzk0OXqwkIRqiDuyS9V9I9oQsCSLrZ3U9197MUDfL8QugCNQHBCJXn7ve7+09DlwNw993uvq7rofskUTPKAcEIAEb3UUl3hC5EE1RlbjoAqBUzu1rSmyW9M3RZmoBgBABDMrM/lnSRpHe5+87Q5WkCghEADMHMrpL0AUWBaHvo8jQFMzCg8uL1qz4k6QRJ/yBpt6QJd/9F0IKhdczsREVrqD2u6FyUpL3M8D0+ghFyF/+H/RNJE5LOlDRL0knu/mTCtm9QtEbVuxWtUfUtRWtUPV1agdFonI/1QDYdirBA0iWStkm6N20jMzta0npJp0i6TNKlkt4k6W4zO6aEcqIdOB9rgD4jFOEed3+tJJnZFZLOS9luuaIxGgvd/Sfx9g9JekzSlZJuLKGsaD7OxxqgZoTcufv+jJteKOm+zn/8+LVPSNoo6X1FlA3tw/lYDwQjhHS6pB8kPL5V0mkllwXgfAyIYISQjlPUjt/rBUnHllwWgPMxIIIRQktK57TSSwFEOB8DIRghpG2Kfo32OlbJv1CBInE+BkQwQkhbFbXT9zpN0g9LLgvA+RgQwQghfUPSud2Lk5nZfEmL4+eAMnE+BsQMDCiEmV0c//Odiqby+bCk5yQ95+4b4m2OkfSgpF2SrlHUXv8JSf9I0lvd/aWyy41m4nysPoIRCmFmaSfWBnd/e9d28zR9+pW7FE2/8mTRZUR7cD5WH8EIABAcfUYAgOAIRgCA4AhGAIDgCEYAgOAIRgCA4AhGAIDgCEYAgOAIRgCA4AhGAIDgCEZARZjZl8zsl/EcaVlf86SZPVlgsUZmZueYmZvZvwldFlQfwQjIyMyuii+u/6qAfU9I+gNJq9395bz3H4K7b5a0VtL1Zvaq0OVBtRGMgOzOjv9uLmDffy7pRUk3F7DvkD4l6XWS/l3ogqDaCEZAdudIeknSY3nu1MzeLOldkr7m7rvy3Hdo7v49SY9IutLMZoQuD6qLYAQMYGZ/ES9BcIqkV0naFzfXuZn9QQ6H+NeKliu4NeX4ZmYfMbOtZrbbzKbM7D+Z2ew+Zf6gma0xs8fNbJeZvWhmG3vLa2anxO9jfZ99PWxme8zsdV2PXWhmd5nZz83sFTN71sw2mNmHE3bxd5LmKQq4QKLDQxcAqIHNkv67pMsk/b2k/9P13IYc9v8uSfsk3Zfy/F8qaub6uaRbJO2R9D5JvynpCEm/TnjNzYqWyr4nft2rJb1X0pfNbKG7/wdJcvdHzOxuSe8wsze7+4+7d2Jmb5N0hqQ17v6L+LEVkj4v6ReSJiU9L+k1kt4q6XJJN/WUZWP8992S1g36MNBS7s6NG7cBN0krFK38uSLn/R4jaa+kh1Oef1t83J9IOq7r8aMkfTd+7smE170x4bEjFC0Wt0fS3K7HL4738+mE13wxfu7dXY9tlvSKpNckbH98wmOz4318L/T3yK26N5rpgGw6yQvfz3m/cyXNUFR7SXJ5/PeT7v5C50F33y3p6rSduvtPEx77taS/VtQi8s6up9ZKelbSB83syM6DZjZH0iWSfirpWz2726soqPUe4/mEx3ZI2q2oqQ5IRDACsjlb0cX34d4nzGx9Vx9S9+3ODPt9dfx3W5/jSsnNgfcqCgqHMLN5ZvbXZvaIme3slEnSmniTuZ1t3X2vpC/EZbmoazeXSpol6RZ3714S+iuSjpa01cw+a2bLzOyE9LcoSXpB0vEDtkGL0WcEDGBmh0t6i6QfuvsrCZsskvRxSf+15/GdGXbfyZ47KuX5TpLC/+t9wt33mdmvEsp7sqTvSTpWUcD635J2KOqXmq+o7+vInpfdIuljkq6U9D/ix1Yo6o/6bz3HvdHMnpf0YUV9WX8oyc1sg6SV7r4p4X3M6nqvwCEIRsBgpykKFoc00ZnZGyXNkbTB4w7+If0y/vvqlOd3xH9fK+nxnmPPiF831fOaq+LHL3f3L/a85gOKgtE07j5lZpOSfsfMTlUUyM6QdKu7P5ew/ZckfSluynubpN9RlBW4zsxOdffO+5KZHaboM3oi5T0CNNMBGZwV/92S8Nw5imocSc9l8XNJz0lamPJ8JwD+y4Tn/rmSf1AuiP+uSXguaT8dnSy4FfFNirLmUrn7dnf/prsvV5TscFxcrm4LFaWuP9BvX2g3ghEwWKfW8mLCc+coSkD4pZm91HW7LcuO476YeyQdb2YLEjb5Yvz342Z2XOdBMztK0ewGSZ6M/769+0EzO1/SFX2Kc5ekHyuqOV0i6cfufnfvRmb2W3HTZa/XxH97myfPjf8esi+gg2Y6YLDO9D+fNLMzJL0saau736YoGK2RtKrnNTuU3RpFiQPnK0rhPsDdN5rZ5yR9VNIPzOx5SZp8AAABx0lEQVR2HRxntE3JWXg3KcrCu83M1ihqxjtD0m9J+pqk30sqhLu7mf1nSTfGD6XViv5O0m4z+46iwGeKakP/RNFn1Zt5d56i2uMdKfsDZNOTZAAkMbOPKAoIv6Go8//P3f3jcQLB9e7+2TH2fYSkpyU95e6/mfC8Sfq38e1kSb+S9D8VJRw8KEnuPr/nNW+TdL2i5IrD4+0+LWm7ohrKde5+bcKxjlU0iLUzFikpQeJDigLnmYrmndst6SlJX5V0s7v/Q9e2sxUNjl3n7ssyfiRoIYIRMCIzO0lRUsE73P3bY+7rakWTpZ7t7qP2P43NzN6uKFj9rbtfmsP+PirpryT9C3e/d9z9obkIRsCIzOxiSbdJOlVRjaPb8/H4naz7OkrSo5Iecvel+ZVyOGb2TUnvkXSuu98/5r5mKRow+/fufnEe5UNz0WcEjO6c+O+Peh53RVllvQEqlbvvNrNLFc0Rd4yXuKaRmb1F0m8rej/vkfS/xg1EsfmKxi99MYd9oeGoGQEtZ2YfVDSw9UVFE5l+OGlaH6BIBCMAQHCMMwIABEcwAgAERzACAARHMAIABEcwAgAERzACAARHMAIABEcwAgAE9/8B2krzU6OAq1oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(t['t_E'], t['pi_E'])\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel('$t_E$ (days)')\n",
    "plt.ylabel('$\\pi_E$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {}
   },
   "source": [
    "## dL vs dS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'$d_S$ (kpc)')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAErCAYAAAAWmx4+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+8HHV97/HXmxDhoNWDgEqOhkDVqJRK4FS9BS1gJQIWI+KvWqlaRNGr1tumBu21/iY19fJDRKSlRVtU5FeotTaoQdAorSdESqlEsTTIASUCB6QcJITP/WN2Yc9mZ3/Mzu7MnH0/H4/zmGR2duYzm5P57Pe3IgIzM7M87FR0AGZmNn84qZiZWW6cVMzMLDdOKmZmlhsnFTMzy42TipmZ5cZJxczMcuOkYmZmuXFSMTOz3DipmJlZbpxUzMwsN04qZmaWGycVMzPLjZOKmZnlZueiAyiLPffcM5YsWVJ0GGZmhdu4ceMvImKvLO91UqlZsmQJU1NTRYdhZlY4SVuyvtfVX2ZmlpvSJRVJT5X0KUnfk3S/pJC0JOXYZ0u6SNIvJM1K2izp3cON2MzM6kqXVICnA68G7ga+nXaQpEngX4FdgBOBo4FPAguGEKOZmbVQxjaVqyPiyQCSTgSObD5A0k7A54BvRsQrGl66cjghmplZK6VLKhHxcBeHHQY8B3jbYKMxM7NelLH6qxuH1ra7SrpG0jZJd0g6U9JYoZGZmY2w0pVUurSotr0QOAtYBUwCHwaeBrwi5X1mNkRrN02zZt1mbpuZZdH4GCuXL2XFsomiw7IBqmpSqZew/iEiPlD787ckLQBWS3pORPxn4xsknQSclHbCxYsXDyZSsxG1dtM0p1x6PbPbtgMwPTPLKZdeD+DEMo9Vtfrrztr26037r6htD2x+Q0ScGxGTaT977ZVp8KiZpVizbvMjCaVudtt21qzbXFBENgxVTSo31LbRtF+1bTeN/WY2QLfNzPa03+aHqiaVrwG/Al7atH95bev5VswKtmi8dZ+ZtP02P5SyTUXS8bU/HlzbHiVpK7A1Iq6KiDslnQr8X0n3AutJGuo/AHwuIm4aftRm1mjl8qVz2lQAxhYuYOXypQVGZYNWyqQCXNT097Nr26tIxqhA0tPrl8DbgT8FbgfWAB8ZQnxm1kG9Md69v0aLIpqbJUbT5ORkeJZiMzOQtDEiJrO8t6ptKmZmVkJOKmZmlhsnFTMzy42TipmZ5cZJxczMcuOkYmZmuXFSMTOz3DipmJlZbpxUzMwsN04qZmaWGycVMzPLjZOKmZnlxknFzMxyU9ap783MKmPtpmlP8V/jpGJm1oe1m6bnLEY2PTPLKZdeDzCSicXVX2ZmfVizbvOc1S0BZrdtZ826zQVFVCwnFTOzPtw2M9vT/vnOScXMrA+Lxsd62j/fOamYmfVh5fKljC1cMGff2MIFrFy+tKCIiuWGejOzPtQb4937K+GkYmbWpxXLJkY2iTRz9ZeZmeXGScXMzHLjpGJmZrlxUjEzs9w4qZiZWW6cVMzMLDdOKmZmlpvSJRVJT5X0KUnfk3S/pJC0pMN7Tqkd953hRGlmZq2ULqkATwdeDdwNfLvTwZL2A94P3DHguMzMrIMyJpWrI+LJEXE0cFEXx38GuAD44WDDMjOzTkqXVCLi4W6PlfT7wEHAKYOLyMzMulW6pNItSbsDpwF/FhF3FR2PmZlVOKkAa4AfAecXHIeZmdVUcpZiSS8ETgAOiojo8j0nASelvb548eKcojMzG12VTCrAZ4HzgFsljdf27QwsqP19NiJ+1fiGiDgXODfthJOTk10lJzMzS1fVpPLs2s/bWrx2N/Ae4PShRmRmZpVNKoe32Hc6sAB4J3DTcMMxMzMoaVKRdHztjwfXtkdJ2gpsjYirIuJbLd4zA+zc6jUzMxuOUiYVdhz0eHZtexVw2HBDMTOzbpUyqUSEMrznsAGEYmZmPajyOBUzMysZJxUzM8uNk4qZmeXGScXMzHLjpGJmZrkpZe8vM7OqW7tpmjXrNnPbzCyLxsdYuXwpK5ZNFB3WwDmpmJnlbO2maU659Hpmt20HYHpmllMuvR6gbWKZD4nI1V9mZjlbs27zIwmlbnbbdtas25z6nnoimp6ZJXg0Ea3dND3gaPPlpGJmlrPbZmZ72g/ZElEZOamYmeVs0fhYT/shWyIqIycVM7OcrVy+lLGFC+bsG1u4gJXLl6a+J0siKiMnFTOznK1YNsGpxx3AxPgYAibGxzj1uAPaNrpnSURl5N5fZmYDsGLZRE89t+rHVr33l5OKmZXGfOhS249eE1EZOamYWSlkHdth5eI2FTMrhfnSpXbUOamYWSnMly61o85JxcxKYb50qR11PSUVJV4i6ROSrpF0m6QHJd0j6ceSLpJ0siRXgJpZT+ZLl9pR11VDvaTdgHcBbwUWA/U15B8A7gDGgP2AXwdeCZwh6SvAJyPiu3kHbWbzz3zpUjvqFBHtD5DeBHwU2Bu4EfgSsAH4fkTc23CcgKXAC4DlwMuBXYCLgZURccsgbiAvk5OTMTU1VXQYZmaFk7QxIiazvLebksp5wFrg1Ij4ftpBkWSnG2s/50t6PPCHwCrgjcCHswRoZmbV0U1SmYyIa3s9ca0U8ylJfw0s6fX9ZmZWPR2TSpaE0vT+B0hKL2ZmNs+5S7GZmeWm56Qi6VWS1ktalPL6hKRvSjqu//DMzKxKspRUTgTGI+K2Vi9GxDTw+NpxZmY2QrIklQOATn1vp4DfzHBuMzOrsCxJ5YkkAx7buRPYM8O5kfRUSZ+S9D1J90sKSUuajpmUdK6kG2vH3CLpAkn7ZrmmmZnlI0tS+QXwjA7HPAOYyXBugKcDrwbuBr6dcsxrgf2BM4GjSMbCHARMSXpaxuuamVmfsqynsgE4VtKzImKHrsKSnk0ymv4rGWO6OiKeXDvXicCRLY75y4jY2nTdDcDNwFuAD2S8tpmZ9SFLSeWvSJLRdyS9S9IzJT22tn03SeliQe24nkXEw10cs7XFvi3AVsATBZmZFaTnkkpEfF/S24FPA6fVfhptB06OiH/NIb6u1UpITwJ+OMzrmpnZozItJxwRfy3pO8DbgecD4yRtKNcAn4mIoT7YJe0MnENSUjlvmNc2M7NHZV6jvpY43pljLP04C/ht4JiIuLvVAZJOAk5KO8HixYsHFJqZ2ejInFTKQtKpJMniDyPiirTjIuJc4Ny01ycnJ9uvAWBmZh1lTiqSHge8AlgGPAG4B9gEXBYR9+UTXscY3k/SnfhdEfH3w7immZmly5RUJL2KpA1jnEdXgQQI4HRJb42Ii3OIr10M7yJZPOz9EfGpQV7LzMy603NSkfQS4IvAw8DngW8BPwOeAhwO/D7wRUkzEfGNLEFJOr72x4Nr26MkbQW2RsRVkl4LnA78C7Be0gsa3n5vRPxnluuamVl/Oi4nvMMbpG+TjF5/Yau1ViRNAlcDUxHxokxBSWlBXRURh0k6n2RVydRjer2mlxM2M0sMejnhZsuAC9MW74qIKUlfBo5v9Xo3IkIdXn8jyRLFZmZWIllG1P8KuL3DMbfVjjMzsxGSJal8Gzi0wzGHkFSBmZnZCMmSVN4LHCBptaTHNr5QmwPsE8BvkHT1NTOzEZKlTeW9wL8DK4GTJF0L/Bx4MkkD/hNISinvleY0jURE/FF/4ZrZqFu7aZo16zZz28wsi8bHWLl8KSuWVWse2flwD2my9P7qOItwioiIBRnfO3Du/WVWfms3TXPKpdczu237I/vGFi7g1OMOqMxDuQr3MOzeX15d0cwKsWbd5jkPY4DZbdtZs25zaR7IncyHe2gnS1K5PSIe7HSQpCUR8d8Zzm9m1tJtM7M97S+j+XAP7WRpqP9CpwNqS/quz3BuM7NUi8bHetpfRoO6h7Wbpjlk9Xr2XfVVDlm9nrWbpvs6X1ZZkspxks5Ie1HS3iQJpfrlODMrlZXLlzK2cG7T7NjCBaxcvrSgiHo3iHuot9NMz8wSwPTMLKdcen0hiSVLUjkL+N+S/rT5BUlPIkko+wIn9BmbmdkcK5ZNcOpxBzAxPoaAifGxUjVwd2MQ99CunWbYsrSpvBtYBKyWdGtEfAlA0h7AN4BnAm+KiAvzC9PMLLFi2USlkkgred9Dmdppei6pRNIH+fXA94DzJR0m6QnAFSSDHk+OiM/nG6aZmaUpU1tTluovIuJXwMuBm4HLSKq8lgHvqa2waGYDUJbGWCuXMrU19bNG/V2SlpOUWA4ETomI1AZ8M+tP86C5emMsUPnqIOtP/d+/DKP0OyYVSX/b4ZAtwGOApU3HeloWsxxlHTQ3n6cEyVPVP6eytDV1U1J5Y5fnaj4uACcVs5xkaYx16aY7/pzy001S8bQsZiWwaHyM6RYJpF1j7HyaEmSQJYn59DkVrWNSiYgtwwjEzNpbuXxpy4kI2zXGDqur6aCrjgZdkihTl9yqy9T7y8yGL8uguWF0NR3GaO5BD+4rU5fcquumoX4sIvpK13mcw8x6b4zNUrrpVTdVR/2WZPIsSbSKZRif06jopqRys6R3S9ql15NLeq6ky4EdpnQxs+HYZedH/5vvvtvC3Kc16fTAb1WSec+FP+DP117f9TXyKkmklaqAyk//UhbdNNRfAfw/4C8kXQh8GbgmreQhaT9gOcncX88DfgqsySdcM2uWVgpotRjUA9uyrrGXrlMHglYlmQAuuOYWJvd5YlcP7rxKEu1KVRtWHeEkkoOOJZWIOAF4PjAFnEQyv9c9kq6T9C+SvijpMklXS/o58GPg08AS4P3A0oj4zsDuwGyEtWvPGNYkg51Gc6eVZAK6jiWvSRjdID94XY2oj4gp4EhJzyAZe/JiklH0BzQduhW4FLgEuCQituUYq5k1aZc4hvUA7TSaO60k02sseQzuy9It23rT0zQtEfFjYBWApN1I1kzZA5gF7oiI23OP0MxSpT2U0x7iMJgHaLsH/srlS3nPhT8ghhRLO26QH7zMXYoj4v6I+HFEXBMR1zmhmA1frw/lIh6gK5ZN8PoXLEYlicUN8oOlZCZ7m5ycjKmpqaLDMOtJq8b4NBMFz2dV9bm1RomkjRExmeW9mWcpNrPiNbdnpH1FFLBh1RFDi6uVskx4aINVuhH1kp4q6VOSvifpfkkhaUmL43aVtEbS7ZJma8e/aPgRmxVrxbIJNqw6gptXH8OER4ZbwUqXVICnA68G7ga+3ea484C3AB8AXgbcDqyTdODAIzQrqTIt1mSjqYzVX1dHxJMBJJ0IHNl8gKTnAr8PvDki/q627yrgBuDDwLHDC9esPMq0WJONptIllYjoZsjvscA24MKG9z0k6UvAKkm71JY8Nhs53bRdVLnRvOjYi75+2ZUuqXRpf+DmiLi/af8NJKtQPr32ZzNrUuUFqXqJfRAP/0F9dnnGWnTS67tNRdJvFdCO8USSNpdmdzW8bmYtDGv6lmZrN01zyOr17Lvqqxyyen2mqfG7jX1Q0/EP4rPLM9ZhLEPQSR4N9auB1zbukPQGSesk/b2k/XO4RjNBy96TzeOrGmM6SdJU2s/WrVsHEKZZ/vp9OBcx/1VeD7tuYx9U4hzEZ5dnrEV9YWiUR1L5DeDy+l9qjeh/R7IM8e8AG1p1Ce7TXbQujeze8PocEXFuREym/ey11145h2iWvzwezkUsSJXXw67b2AeVOAfx2eUZaxkmzMwjqTweaPyN/gPgRmApsB+wgdp8YTm6Adi3Nv9Yo+cADwI35Xw9s1LI4+HcqtsxwP0PPjSwapK8Hnbddpke321hy/fvJPV1j4Posp1noirDCpZ5JJWfkkwsWXcEcHEkHgI+ARyew3Ua/SOwEHhVfYeknYHXAFe455fNV3k8nOvzX42PzX3w3n3/toHVv+f1sOtm7q61m6a574GHWr5/e0Rf9ziIucPyTFRlGKeUR++vK4CVwHG1BbqeC7yr4fWbgaf1ckJJx9f+eHBte5SkrcDWiLgqIn5QWzDsdEkLa9c4maTK7fXZb8Ws3PKaun3FsgnWrNvMzOzc1SmalwHOS56zA3fqMr1m3Wa2PZw+p2G/95j3dDN5ji0qwzilPJLKx4FNkqZJuvNuAb7b8PrewC97POdFTX8/u7a9Cjis9uc3AR8DPgqMA9cBL42Ia3u8llll5PlwHmb9e6eHXZ7dYLuJv2yLcuWZqIqeY63vpBIRt0n6LeDdJA/3M2Pu1McvBn7U4zlTe3E1HDML/J/aj9lIyPOb6LAXrEp72OU99uMJYwt3KIE16+Uee0l4RY8RKYOOSUXSxcBG4Frg2ojYoe9tRNwC/EnKKZ4NXNxPkGb2qLy+iZZlwap2nQ+y3Kc6fCXt5R57HWxZ1UGleeqmpHJc7ScAatVc19Z+NpIkmtQFuiLiDTnEaWY5K0P9O+RfDTdzf/tSyi47d98/qZeEl3dyrKpukso+wEG1n4Nr22OB36sfIOnn7Jhofpp7tGaWq6Lr3yH/ari089XNzG7rugTRS8IrwxiRMuiYsiPipxFxeUT8RUS8LCIWAX9DMnr9RpJxKA8DRwPvBy4F/ruWaMyshPKYMiUveXeDXbl8KQt3al8H1u3Ynl66QpdhjEgZ9DxORdJKkm67L4qI/SPiRRHxVOBQYD1JsrkVGK30bFYRg54fqteElffYjxXLJnjcrp0rYbopQfSS8MowRqQMsvT+egfwxYj4TuPOiPgu8BJJ7yVZOOuQHOIzy2xUeuL0ep+DrPvP2liddzVcp3YV6K4E0Uu7U1naqIqWJak8Cbgz7cWI+EtJxwHvI1nB0WzoRqUnTpb7HGTdf1kaqzu1q/RSgugl4ZWhjapoWaZp+THwux2O+RbwwgznNstFGWZrHYYs9znIuv+yNFa3qoqqt7LkMbWKpctSUvlb4DRJ74uIj6cc8xS8pokVqCwPt0HLcp+DHJ8y7AGVaYZZFTUq1azdypJUzgKOAT4i6XeBDzS2r0g6hmR9FXcptsKU5eE2aFnuc5AP3LIMqIThVEWNSjVrL3pOKhGxvZY4zgb+CLhK0p3ALSTtLRMkJc0z8gzUrBdlergNUtb7HNQDt9eEVfVldMvShlQmmjtNV49vlp4HvBN4CUlC2U7S5nJqRPx9LhEOyeTkZExNTRUdhuWo+SFz+LP24sobt867aooqVL+0ihFomRCztHc0lxj6OVe38a9YNsG+q76augTtzauP6fu6RZG0MSImM723n6TSFMQuwPbaGiqV46Qyvw36oVMGZU0uaZ/9rgt34u4WXX8nxsfYsOqInq5xyOr1qb29JnIoAaX97qxZt7nldbPcQ5n0k1TymPoeAC+MZWU236spyli3X09yrR66s9u27/DvUZfnMrrQ+2fRnJzvf/Ch1N+dUalm7UVuScWszKrcG6ybEkgeSTPv9o3mh223si6j225cSmM363b32Co5p7ltZtYDHltwUrGRUNXeYN2WQPpNmnmXdFoluWZjC3figW0Pz2mT6GcZ3U5JbHpmlj+56Dq211aFnJ6ZZeXF1wFzOxh0mwjrvzse8DhXHmvUm5VeVedl6nZwY7cDGtPm5cp7sGinZLZwJ/HQwzEnoQh45cHZHtCN84e1s71pmeFt24MPfeWGR/7ebRKuwu9OUZxUbCTkPWnhsHRbAukmababSLLb63Q7WWS7EuDE+BiP23Vntm2f+4AP4Mobd1gDsGsrlk2wYdURnP6aA3f4LNpp7CyQFvf42MLK/e4UxdVfNjKqWE3RbbVdN3X77Uoj3VynlyqytAbs+sN431VfbXm/ebRx1WP54wt/0PN70+L+4LH7V+53pyguqZiVWC/VdvVv6jevPoYNq47Y4SHYrjTSzXV6qSLrVDIc9NojK5ZNdKwKqxsfWzjnfVUs0ZaJSypmJZZn76K00sj4bgvnXGd6ZpYF0pyEsWLZRM+dAdqVDFuVCAQc/qy9ermltrppvF+4k/jgsfvP2TeoEm1ZxxHlzUnFcjFf/8OU4b7yesitXL6UlRdft0Nbxn0PPMTaTdOPXCOtiiuvHnT1z7T5YR/AJRunmdznibncb6uEXNSsCmUcRzQouY2orzqPqM9uvo5Wn4/3deCHrmBmdsdR7AskPvnq57YdId6pnaQb3Yxf6XY0ehkSfrfSRvyXdeR9PyPq3aZifZuva5f0e19lWge+rlVCAdge8UjPsFbqA/36bW/oZhxIN431g14SOW9VHnzbK1d/Wd/m63+Yfu6rjNUdazdNI2g5ASIkCXOBxPYWtRfdDvTrVHro5rNrV53WaeqXsk67U9XBt1m4pGJ9G3RPnl7kWTro577KWHpbs25zakKp2x7xyAqJdd0O9Oum9NDps2t3rcbzpynrF5mqDr7NwknF+laW/zB5V4n0c19lLL11c+3mkkzjKPdOCbubRLpy+dIdklbdAqltdVo3VWdl/eY/Sl2VXf1lfSvLpHp5z0Tcz32Vsbqj06SLrarG6qPc06rzprbc9UhvqrRSUGMyW7Fsgqktd3HBNbfsMOdXp4dsp6RY9m/+VRx8m4WTiuWiDP9hBlE6yHpfZZwSPW1sSJB8c27XSJ+WsJuTQyvNifSjKw5gcp8n9pysx3db2HL9Feh/zRTLT2WTiqRDgL8ADgR2BW4CzoqIvy00MCtMmUoHZSm99RJTWrfXnaTUhNMpobQb/d/rZ5E2+mF8bGEpu+WOqkomFUm/CXwDuAZ4C3A/cDxwnqRdIuIzRcZnxShb6aAMpbdmvY5yB1r2ButEkHsivSelO3TafitGJZMK8FpgAfB7EXFfbd/XJT0XOAFwUhlBZSwdVEnz57dTSvfiurTuyYMa0JdWEt1JmjMjgBWrqknlMcA2oPk3bAbYffjhWFmUsXRQJY2fX9pMwpAkjsOftReXbJweWsmwXUmq6DFA9qiqdik+v7Y9U9IiSeOS3gK8GDituLDM5o+0tqh6SeSjKw4YajfZerfcBdqxU3LRY4DsUZWd+0vSbwGXAfXf4G3AyRFxXsrxJwEnpZ1v8eLFB2/ZsiX3OM3yMuy5rso699m+q77astpNwM2rjxl2OPNSP3N/VbL6S9IzgEuAG4C3kVSDvRw4R9IDEXFB83si4lzg3LRzTk5OVjO72kgoYtqXQbdRZU2SZerlZzuqZFIBPk5SMnlZRNS7fnxT0h7AGZK+GBEPFxeeWb7yHtjZraxtVJ0SRj9Jsmy9/GyuqrapHABc15BQ6v4N2AN40vBDMhucMk77kqab6XL6mRttlKY8qaKqllR+Bhwo6TER8WDD/ucDDwB3FROW2WCUtcqnVYmkm1JVv0nSvfzKq6pJ5SzgIuArks4maVM5FngdcFpTojGrlFYP6jJW+aRVYaVN+tiYMMqaJK1/laz+ioiLgaOBXYC/IWm0PxR4B7CywNDM+pJWdQSUrsonrUTSqssvzE0YRcxsXcZF0+ajqpZUiIivAV8rOg6zPLWrOtqw6ohSVfmkVVVtj2Bs4YK2paphz35QxkXT5qvKJhWz+WjQDfJ5jnVJq8KaaGhbaXedYbaLFNV7bhQ5qZiVyCDbGvL+tt6unadsDelV6j1XdZVsUzGbrwbZ1tBNN95e2h2q1LW3TEtez3cuqZiVyCDbGjp9W89SkilbiSRNGXvPzVdOKn0Y9lxMNhoG9aDuVLU2n9sdvCzC8DipZOTeJFY1nb6tp5VkpmdmOWT1+so/hKtSqqo6t6lk1M80E2ZF6NQG0q59odVUK2atuKSSkXuTWBVlWU64Lo+qMFcZz39OKhl5mgmbbxrbHVr9bkN/X5pcZTwaXP2VURHTTJgN2oplE2xYdQQTGbvgtuuS3GuVsadVqSaXVDJybxIrwrCqj7J0wW1XEqn/vZVWpR+XaqrLSaUP7k1iwzTMB22WL01pJZEP/uMN/Oqh9DXzWpV+5nP35vnOScWsIob9oO31S1Nae8vMbPNaeo9KK/24I0x1uU3FrCLK/qDN0kklbVoXT6tSXU4qZhVR9gdtWueV3Xdb2PL4ifGxtt2b3RGmmlz9ZVYRw5y/KkuHgLR2GKDnuN0RproUEUXHUAqTk5MxNTVVdBhmbQ2j91dzhwBIkkA/MxB70GO1SNoYEZOZ3uukknBSMUscsnp9X4tvWfX1k1Rc/WXWZNS/VbebWNJjR6wTN9SbNahX/UzPzBKM5kSKaQ3/CyRPomodOamYNfDs0+k9r7anVJWXpUuzlYOTilmDso8FGYa0KfKzzgdmo8VtKmYNPPt0Im00vZfktU5cUjFr4EF36Tot8mUGLqmYzeFBd+15ElXrxEnFrIkfnGbZufrLzMxy45KK2QCM+gBKG12VL6lIOlrS1ZLuk3SvpClJRxQdl40uD6C0UVbppCLprcDlwEbgFcCrgIuA3YqMy0bbIAdQet12K7vKVn9JWgKcDqyMiNMbXlpXSEBmNYMaQOl1260KqlxSeTPwMHBO0YGYNRrUYlqeQsaqoMpJ5VDgRuC1kn4i6SFJN0l6R9GB2Wgb1ABKTyFjVVDZ6i9gUe1nDfA+4CckbSpnSdo5Is4oMjgbXYMaQOkpZKwKKrtIl6QfAc8AXhkRlzbs/xqwDNg7Gm5O0knASWnnW7x48cFbtmwZYMRm/RnEioxmrYzqIl13kiSVrzftvwJ4KbA3cFt9Z0ScC5ybdrLJyclqZlcbGZ5CxqqgyknlBuAFLfartn14iLGYDYWnkLGyq3JD/WW17fKm/cuBWyPiZ0OOx8xs5FW5pPLPwJXAZyXtCfwXcDxwJPCmIgMzMxtVlU0qERGSVgCnAh8CdifpYvz6iPhCocGZmY2oyiYVgIi4F3hH7cfMzApW5TYVMzMrmcqOU8mbpK3AfB+osifwi6KDKCl/Nu3580k3Hz+bfSJiryxvdFIZIZKmsg5omu/82bTnzyedP5u5XP1lZma5cVIxM7PcOKmYmVlunFTMzCw3TipmZpYbJxUzM8uNk8poSZ363/zZdODPJ50/mwYep2JmZrlxScXMzHLjpGJmZrlxUpnHJB0v6RJJWyTNStos6VRJv1Z0bGUk6V8khaSPFh1LWUg6WtLVku6TdK+kKUlHFB1X0SQdIukKSXfUPpdrJb256LjKwEllfvtTYDvwPuClwGeAk4GvS/K/fQNJrwOeW3QcZSLprcDlwEbgFcCrgIuA3YqMq2iSfhOFqVoAAAAIIUlEQVT4BrAQeAvwSuD7wHmSTi4ytjJwQ/08JmmviNjatO8E4HPAiyNifTGRlYukcZIF3t4DfAH4WET8ebFRFUvSEuCHwCkRcXqx0ZSLpI+TfGF7YkTc17D/GpL1A/9XYcGVgL+tzmPNCaXm+7XtxDBjKblPADdExBeLDqRE3gw8DJxTdCAl9BhgGzDbtH8GP1P9AYyg36ltf1hoFCUh6VDgBODtRcdSMoeSlN5eK+knkh6SdJMkr7IK59e2Z0paJGlc0luAFwOnFRdWOVR6OWHrjaQJ4MPANyJiquh4iiZpIfBZ4K8iYnPR8ZTMotrPGpI2uZ+QtKmcJWnniDijyOCKFBH/Iekw4DIe/TKyDXhbRHypsMBKwkllREh6HEmj60PAmwoOpyzeC4wBHys6kBLaCfg14I0RcWlt3/paW8spks6MEW2QlfQM4BLgBuBtJNVgLwfOkfRARFxQZHxFc1IZAZJ2Bf4R2A/4nYi4teCQCidpMfB+4ERgF0m7NLy8S63x/pcRsb2QAIt3J/AM4OtN+68g6Um4N3DbsIMqiY+TlExeFhHbavu+KWkP4AxJX4yIh4sLr1huU5nnalU8lwDPA46OiOsLDqks9gN2Bf4BuLvhB5KePXcDBxQTWinckLJfte3IPjRJfi+ua0godf8G7AE8afghlYeTyjxWG4tyAUkD4ssj4pqCQyqTHwCHt/iBJNEcDtxUTGilcFltu7xp/3Lg1oj42ZDjKZOfAQdKekzT/ucDDwB3DT+k8nD11/z2aZLG1Y8B/yPpBQ2v3TrK1WARMQN8q3m/JIAtEbHDayPmn4Ergc9K2hP4L+B44EjcJncWySDQr0g6m6RN5VjgdcBpEfFgkcEVzYMf5zFJ/w3sk/LyhyLig8OLphokBR78CICkxwOnkiST3Um6GK+OiC8UGlgJSDqKpKPH/iTVqD8hmQL/syPcDgc4qZiZWY7cpmJmZrlxUjEzs9w4qZiZWW6cVMzMLDdOKmZmlhsnFTMzy42TipmZ5cZJxczMcuOkYtaGpMWSQtKlnY/u+pyfl3SHpMc27FtSu875eV0nK0kH12L5o6JjsepxUjFr76Da9to8TiZpEvgDkulO/iePc+YtIjYCa4GP1tbhMeuak4pZe/WksjGn830cuBf4TE7nG5RTgacA7yo6EKsWJxWz9g6ubfsuqUh6JvC7wJcjYrbf8w1SRPwbyQSSb5W0oOh4rDqcVGzkSdpZ0rsl/bukWUlbJP2ZknnwDwKmI+LnOVzqzSSLXF3YQ2w7STqz3q4jadfG9hdJz5K0VtJdkv5H0nckHdnmfM+TdKGkaUm/knS7pCskvbrF4V8CFpMkQrOuOKnYSKsttPTPwOnAdpK1Mr4JfIhkKvOnkFN7CsnDeTvQ1WJptWWgvwy8k2RtnOMj4oGGQ/YFvkey2uBnSdb4OBj4mqTXtDjfW4DvAitq208CXyVZqfDtLULYUNu+pJt4zcCLdJl9muSh+QHgo1FbC6LWC+uq2jF5VH09FjgQ+GE3DfSSnghcDhwCrIqIv2xx2IuAv4qIlQ3vO4sk0Zwj6WsRcW9t/3OAs0nac14YEXOWC5b01Bbn/37Ddcy64pKKjSxJzwNOBK6OiI9Ew+JCEXE18MPaX/NopJ8AFgC3dxHXPiSlhOcDb0hJKAD3AB9u3BERUyRLSI8Dr2h46WSSL5EfaU4otfftsApoRNxDsjzu4k4xm9U5qdgoe2dt+4GU1++sbR8pqUhaL+nzGa61R217d4fjlpKUNBYBR0XEBW2OvTYiftli/7dq22UN++pLSX+tw/Wb3QXs2eN7bIQ5qdgoO5IkcVyd8vp+wM8jYrph3zKylVzqvb127XDcM4G9SdaE71TtltZ54Ge17RMa9o3XttP0ZoxHYzfryEnFRlKtEfxJwC3RYk1tSb9NUlpoLKX8OsnDOUtSuaO23aPtUfAV4H0k7S/flNSulPDklP1PqW3vadg3U9tOdLj+IyTtRHK/d3Q61qzOScVG1fbaz5NSXv9QbdtYWjgYeBj4QYbr3Q5sJaneaisiTgXeQ1IqulJSWvI4SNKvtdh/WG27qWFfvcfZUV1Fm1hK0gU6y/3aiHJSsZEUEduAHwMTkn6v8TVJ7+XRsRmNpZKDgR9FxH0Zrhck1Wx7Snp6F8efTtK4vj9wlaRFLQ57Ak3tQbVpYF5PUkq5rOGlzwAPAf+31hOMpve16v1Vb4e5slO8ZnXuUmyj7FTgc8Alkr5E0hZxGHAA8FPgaexYUumnJ9glwCuB5cBNnQ6OiHMkPQCcB1wt6YiIuKXhkKuBEyU9n6S32N7Aa0i+LL613p24dq7/lPR24Bxgk6TLSZLqHsAk8Evg8KYQjiQpzV2e5WZtNLmkYiMrIj4P/DFwK/A64A9JkskhQAB3RsSWhrdkbaSvu4Skcf2EHmI8n2QCyn1IEst+DS/fDPw2SY+ytwGvJkmCR0fEDqP2I+KvgUOBfyJJniuBY4FfkIzXeYSkJ5AMkvyniPhpt/GaqUUbpZk1kbQvSY+swyLiqk7HtznPKSSTSh4UEZs6HZ9yjiUkCeVzEfHGrLF0uMY7gTOBF0XEtwdxDZufXFIx6059Ysltkn6j4efZPZ7nNOAWmgYtlomkMeAU4BInFOuV21TMulNPKhua9v8HSRtMVyLiAUlvAA6X9NiSrqmyhGTes/OLDcOqyNVfZhUzjOovs6ycVMzMLDduUzEzs9w4qZiZWW6cVMzMLDdOKmZmlhsnFTMzy42TipmZ5cZJxczMcuOkYmZmufn/UtXHf6kjiWwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(t['rad_L'], t['rad_S'])\n",
    "plt.xlabel('$d_L$ (kpc)')\n",
    "plt.ylabel('$d_S$ (kpc)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {}
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.5"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "source": [],
    "metadata": {
     "collapsed": false
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
